//这是一个自定义内核的图染色算法，主要针对的是HPCG中的线性方程组的系数矩阵A进行图染色。
#include<stdio.h>
#include<stdlib.h>
#include<string.h>
#include<cusparse.h>
#include<sys/time.h>


#define A_row_size 32768

__global__ void init_need_to_color_device(int* need_to_color_device, int size){
    //每个线程为一个位置附上一个数字
    int thread_id = threadIdx.x + blockIdx.x * blockDim.x;

    //附上一个数字
    if(thread_id < size){
        need_to_color_device[thread_id] = thread_id;
    }
}

//我们写一个函数来记录每个颜色的节点数量
__global__ void get_node_num_of_each_color(int* coloring, int row_num, int* color_num_record){
    //获取每一个线程的线程号
    int thread_id = threadIdx.x + blockIdx.x * blockDim.x;

    if(thread_id < row_num){
        int color = coloring[thread_id];

        //这里用一个原子加
        atomicAdd(&color_num_record[color], 1);
    }
}

//传入一个二维数组，每行二维数组元素的个数，二维数组的行数，传入一个数组作为输出，标记了每个节点的颜色，另外还有一个count
//计算已经染色的节点数量，color_num是已经有的颜色数量，每次迭代增加两个颜色。最后需要一个输出数组，计算每个节点的颜色
//pitch是数组每一行对齐之后的大小
//必须传入两个coloring数组，有一个数组是要在本次迭代中被修改的，有一个数组记录的是上次迭代的结果，一个节点是不是要在
//这次迭代和自己的邻居节点比较，是取决于这个邻居节点是不是在上次迭代的时候被染色过

//每次只染色需要染色的点，在数据规模过大的时候，对所有节点进行遍历非常低效，所以我们引入need_to_color数组
__global__ void max_min_graph_coloring(int* A_device, size_t pitch,int* element_num_each_row_device, 
    int need_to_color_num, int* need_to_color_num_new, int color_num, int* coloring, int* old_coloring_device, 
    int* need_to_color, int* need_to_color_new){
    //计算当前节点编号
    int thread_id = threadIdx.x + blockIdx.x * blockDim.x;
    
    //如果这个线程负责的行号存在才执行下列逻辑
    if(thread_id < need_to_color_num){

        //通过当前线程编号获取需要染色的节点编号
        int target_point = need_to_color[thread_id];

        if(old_coloring_device[target_point] == 0){
            // printf("开始执行染色\n");
            //存储当前节点邻居的数量
            int col_num = element_num_each_row_device[target_point];

            //计算当前节点的哈希值，哈希函数是 thread_id % 29
            float random_num = target_point % 23;

            //查看所有其他节点的哈希值
            float random_neighbor = 0;

            //用两个值判断这个是极小值还是极大值，假设一开始即是极大值又是极小值
            int max = 1;
            int min = 1;
            

            //遍历所有的邻居节点
            int i;

            int neighbor_col;

            //和所有相邻节点相比较，极大值染色color + 1，极小值染色color + 2
            int* A_device_row = (int*)((char*)A_device + target_point * pitch); 
            for(i = 0; i < col_num; i++){
                neighbor_col = A_device_row[i];

                // if(thread_id == 10){
                //     printf("%d ", neighbor_col);
                // }

                //查看邻居节点有没有染色，如果已经染色了，就不要和邻居节点比
                if(old_coloring_device[neighbor_col] == 0){
                    //计算所有节点随机数
                    random_neighbor = neighbor_col % 23;

                    //计算这个随机数和当前节点随机数的比较
                    if(random_neighbor > random_num && max == 1){
                        max = 0;
                    }

                    if(random_neighbor < random_num && min == 1){
                        min = 0;
                    }

                    //如果两个节点一样，如果都是极小值，那么节点编号比较小的就是比较小
                    //如果都是极大值，那么节点编号比较小的就是比较大
                    if(random_neighbor == random_num && target_point > neighbor_col){
                        min = 0;
                        max = 0;
                        break;
                    }
                }
            }

            // if(thread_id == 10){
            //     printf("\n");
            // }

            //根据实际情况染色
            //如果一个节点又是极大又是极小说明他没有邻居了，按照极小处理
            if(min ==  1 || max == 1){
                if(max == 1){
                    coloring[target_point] = color_num + 1;
                }
        
                if(min == 1){
                    coloring[target_point] = color_num + 2;
                }

                // atomicAdd(count_device, 1);
            }else{
                //这里说明当前节点没有被染色，就加入need_to_color
                int insert_index = atomicAdd(need_to_color_num_new, 1);
                need_to_color_new[insert_index] = target_point;
            }
        }
    }
}

//检查颜色是不是打印正确的内核函数
//传入矩阵A以及每个节点的染色情况，看看某一个节点的邻居节点的颜色是不是冲突了
__global__ void check_graph_coloring(int* A_device, size_t pitch, int* element_num_each_row_device, 
    int row_num, int* coloring){
    //遍历每一行
    int thread_id = threadIdx.x + blockIdx.x * blockDim.x;
    
    if(thread_id < row_num){
        //将一行取出
        int* A_device_row = (int*)((char*)A_device + thread_id * pitch);

        int col_num = element_num_each_row_device[thread_id];

        int neighbor_col;

        int i;

        //获取当前节点的颜色
        int color = coloring[thread_id];

        if(color == 0){
            printf("出现了错误，有节点没有染色\n");
        }

        for(i = 0; i < col_num; i++){
            //查看当前节点的列号
            neighbor_col = A_device_row[i];

            if(coloring[neighbor_col] != 0 && color != 0 && color == coloring[neighbor_col] && neighbor_col != thread_id){
                //这里说明染色冲突
                printf("出现错了，%d与%d冲突，冲突颜色%d\n", thread_id, neighbor_col, color);
            }
        }
    }
}




int main(){
    //我们创造一个二维数组，这个二维数组一共4096行，27列，每一行就代表了图的节点数量
    //每一列就代表了就是一个点，一行中的不同元素代表了这个节点和哪些行有连接。
    //创造4096个指针
    //把A改成堆区
    int** A = 0;

    cudaMallocHost((void**)&A, A_row_size * sizeof(int *));

    //为每个指针创造27个空间
    //遍历变量
    printf("申请空间\n");

    int i;
    for(i = 0; i < A_row_size; i++){
        cudaMallocHost((void **)&A[i], 27 * sizeof(int));
    }

    //用一个堆区数组来记录每一行元素的数量
    int* element_num_each_row = 0;

    cudaMallocHost((void **)&element_num_each_row, A_row_size * sizeof(int));

    cudaMemset(element_num_each_row, 0, A_row_size * sizeof(int));
    
    //我们将数据写入。
    
    //从文件中读入数据
    FILE* Graph_File;

    //打开文件
    Graph_File = fopen("../data/A_Matrix_col_1.txt", "r");

    //用一个缓冲区来获取读取的内容
    char read_line_buf[1024];

    //用一个变量来遍历A的每一行
    int A_row_index = 0;

    printf("读取文件\n");

    fgets(read_line_buf, 1024, Graph_File);

    while (!feof(Graph_File)) {
        
        //用一个变量来定位A中每一行的每一个元素
        int A_col_index = 0;
        
        //分割字符串，这里是每一个子串
        char* part_char;

        part_char = strtok(read_line_buf, " ");


        while(part_char != NULL && A_col_index < 27){

            //第一个元素
            A[A_row_index][A_col_index] = atoi(part_char);

            // printf("%d ", A[A_row_index][A_col_index]);
            
            A_col_index++;
            element_num_each_row[A_row_index]++;

            part_char = strtok(NULL, " ");
        }

        //先读，才能查到eof
        fgets(read_line_buf, 1024, Graph_File);

        A_row_index++;
    }

    fclose(Graph_File);

    // printf("开始打印\n");
    // //打印A矩阵，看看情况
    
    // for(i = 0; i < 4096; i++){
    //     for(j = 0; j < element_num_each_row[i]; j++){
    //         printf("%d ", A[i][j]);
    //     }

    //     printf("\n");
    // }

    //这里打印一下element_num_each_row
    // for(i = 0; i < 4096; i++){
    //     printf("%d\n", element_num_each_row[i]);
    // }

    //开始进行图染色
    
    //首先在GPU端申请二维数组，使用cudaMallocPitch
    //https://www.cnblogs.com/csyisong/archive/2010/01/10/1643519.html
    
    //设备端A矩阵指针
    int* A_device = 0;
    //每一行对齐之后的大小
    size_t pitch;
    
    printf("申请设备端的A指针\n");
    //申请设备端的A指针
    cudaMallocPitch((void**)&A_device, &pitch, 27 * sizeof(int), A_row_size);

    //将数组A的数据拷贝给设备端，值的一提的是，数组A的每一行在存储上不是连续的，所以我们使用memcpy2D来处理
    //我们所能用的方式就是将A矩阵一行一行地传给设备端
    printf("pitch = %d\n", pitch);

    for(i = 0; i < A_row_size; i++){
        //因为pitch是字节数量，如果想让A_device这种int指针做地址偏移的话，需要将A_device的指针类型变化一下。
        cudaMemcpy((char *)A_device + i * pitch, A[i], 
            element_num_each_row[i] * sizeof(int), cudaMemcpyHostToDevice);
    }
    cudaThreadSynchronize();


    printf("将数组A的数据拷贝给设备端完毕\n");

    //然后将每一行多少个元素的数组拷贝到GPU端
    int* element_num_each_row_device;

    cudaMalloc((void **)&element_num_each_row_device, A_row_size * sizeof(int));

    //将数据拷贝到设备端
    cudaMemcpy((void *)element_num_each_row_device, element_num_each_row, A_row_size * sizeof(int), cudaMemcpyHostToDevice);

    printf("然后将每一行多少个元素的数组拷贝到GPU端\n");

    //创建一个变量来保存需要染色的节点数量
    int need_to_color_num = A_row_size;
    int* need_to_color_num_device;
    int* need_to_color_num_new_device;

    //将这个变量传到显存上
    cudaMalloc((void**)&need_to_color_num_device, sizeof(int));
    cudaMalloc((void**)&need_to_color_num_new_device, sizeof(int));

    cudaMemcpy(need_to_color_num_device, &need_to_color_num, sizeof(int), cudaMemcpyHostToDevice);

    //初始化need_to_color
    int* need_to_color_device;
    int* need_to_color_device_new;

    cudaMalloc((void**)&need_to_color_device, A_row_size * sizeof(int));
    cudaMalloc((void**)&need_to_color_device_new, A_row_size * sizeof(int));

    //初始化need_to_color，一开始所有节点都是需要染色的
    //创建线程结构，线性网格，每个block64个线程
    dim3 block_dimension(384);

    dim3 grid_dimension(need_to_color_num / 384 + 1);

    init_need_to_color_device<<<grid_dimension, block_dimension>>>(need_to_color_device, A_row_size);
    
    //创建一个数组来表达每个节点所染的颜色
    int* coloring = 0;
    
    cudaMallocHost((void **)&coloring, sizeof(int) * A_row_size);

    int* coloring_device = 0;

    cudaMalloc((void **)&coloring_device, sizeof(int) * A_row_size);
    //将整块显存初始化为0
    cudaMemset(coloring_device, 0, sizeof(int) * A_row_size);

    //上次迭代的染色情况
    int* old_coloring_device = 0;

    cudaMalloc((void **)&old_coloring_device, sizeof(int) * A_row_size);

    //创建颜色数量
    int color_num = 0;


    //循环调用内核函数
    //暂时只进行一次迭代
    //计时开始
    struct timeval starttime,endtime;
    
    double timeuse;

    gettimeofday(&starttime,NULL);

    while(need_to_color_num != 0){
        printf("开始运行内核函数，need_to_color_num = %d\n", need_to_color_num);

        //将上次迭代的结果暂存起来
        cudaMemcpy(old_coloring_device, coloring_device, sizeof(int) * A_row_size, cudaMemcpyDeviceToDevice);

        // max_min_graph_coloring(int* A_device, size_t pitch,int* element_num_each_row_device, 
        //     int need_to_color_num, int* need_to_color_num_new, int color_num, int* coloring, int* old_coloring_device, 
        //     int* need_to_color, int* need_to_color_new)
        max_min_graph_coloring<<<need_to_color_num / 384 + 1, block_dimension>>>(A_device, pitch, element_num_each_row_device, need_to_color_num, 
            need_to_color_num_new_device, color_num, coloring_device, old_coloring_device, need_to_color_device, need_to_color_device_new);

        // cudaThreadSynchronize();

        color_num = color_num + 2;
        
        //将need_to_color_num拷贝回来，作为下一轮的需要染色的节点数量
        cudaMemcpy(&need_to_color_num, need_to_color_num_new_device, sizeof(int), cudaMemcpyDeviceToHost);
        //将need_to_color_device_new，拷贝到needs_to_color_device，作为下一轮需要染色的节点统计
        cudaMemcpy(need_to_color_device, need_to_color_device_new, need_to_color_num * sizeof(int), cudaMemcpyDeviceToDevice);

        //初始化下一轮需要染色的节点数量
        cudaMemset(need_to_color_num_new_device, 0, sizeof(int));
    }

    cudaThreadSynchronize();

    //计时结束
    gettimeofday(&endtime,NULL);

    timeuse=(1000000*(endtime.tv_sec-starttime.tv_sec)+endtime.tv_usec-starttime.tv_usec)/1000;
    
    printf("使用的时间:%fms\n", timeuse);
    //查看错误
    //__global__ void check_graph_coloring(int* A_device, size_t pitch, int* element_num_each_row_device, 
    //int row_num, int* coloring)
    check_graph_coloring<<<grid_dimension, block_dimension>>>(A_device, pitch, element_num_each_row_device,
        A_row_size, coloring_device);

    cudaThreadSynchronize();

    //将颜色拷贝到cpu端
    cudaMemcpy(coloring, coloring_device, sizeof(int) * A_row_size, cudaMemcpyDeviceToHost);

    //打印所有的染色结果
    // for(i = 0; i < 4096; i++){
    //     if(coloring[i] != 0){
    //         printf("%d:%d ", i, coloring[i]);
    //     }
    // }

    // printf("\n");

    //这里计算每一个节点的颜色
    //这里检查每个颜色的节点的数量
    int* color_num_record;
    int* color_num_record_device;

    cudaMalloc((void **)&color_num_record_device, (color_num + 1)*sizeof(int));
    cudaMemset(color_num_record_device, 0, (color_num + 1)*sizeof(int));
    
    cudaMallocHost((void **)&color_num_record, (color_num + 1)*sizeof(int));

    //每个颜色的节点数量统计
    //__global__ void get_node_num_of_each_color(int coloring, int row_num, int* color_num_record)
    get_node_num_of_each_color<<<A_row_size / 384 + 1, block_dimension>>>(coloring_device, A_row_size, color_num_record_device);

    //将数据传回
    cudaMemcpy(color_num_record, color_num_record_device, (color_num + 1)*sizeof(int), cudaMemcpyDeviceToHost);

    printf("每个颜色的染色数量\n");

    //打印数据
    for(i = 0 ; i < color_num + 1; i++){
        printf("%d, ", color_num_record[i]);
    }

    printf("\n");
    

    cudaFree(color_num_record_device);
    
    cudaFreeHost(color_num_record);

    
    printf("颜色种类：%d\n", color_num);

    //析构A
    
    for(i = 0 ; i < A_row_size; i++){
        cudaFreeHost(A[i]);
    }



    //析构矩阵A
    cudaFreeHost(A);

    //析构coloring_device
    cudaFree(coloring_device);

    //析构coloring
    cudaFreeHost(coloring);

    //析构element_num_each_row_device
    cudaFree(element_num_each_row_device);

    //析构设备端矩阵A
    cudaFree(A_device);

    //析构need_to_color_num_new_device
    cudaFree(need_to_color_num_new_device);

    //析构need_to_color_num_device
    cudaFree(need_to_color_num_device);

    //析构need_to_color_device
    cudaFree(need_to_color_device);

    //析构need_to_color_device_new
    cudaFree(need_to_color_device_new);

    

    //析构old_coloring_device
    cudaFree(old_coloring_device);
}